Cite this Article
Li Chen, Liu Peng, Li Yongfang, Wang Dunyou. Solvent effects and potential of mean force study of the SN2 reaction of CH3+CN− in water. Chinese Physics B, 2018, 27(3): 033401  
Permissions
Solvent effects and potential of mean force study of the SN2 reaction of CH3+CN− in water
Li Chen, Liu Peng, Li Yongfang, Wang Dunyou ‡
College of Physics and Electronics Shandong Normal University, Jinan 250014, China

 

† Corresponding author. E-mail: dywang@sdnu.edu.cn

Abstract
Abstract

We used a combined quantum mechanics and molecular mechanics (QM/MM) method to investigate the solvent effects and potential of mean force of the CH3F+CN− reaction in water. Comparing to gas phase, the water solution substantially affects the structures of the stationary points along the reaction path. We quantitatively obtained the solvent effects’ contributions to the reaction: 1.7 kcal/mol to the activation barrier and −26.0 kcal/mol to the reaction free energy. The potential mean of force calculated with the density functional theory/MM theory has a barrier height at 19.7 kcal/mol, consistent with the experimental result at 23.0 kcal/mol; the calculated reaction free energy at −43.5 kcal/mol is also consistent with the one estimated based on the gas-phase data at −39.7 kcal/mol.

PACS: 34.20.-b;34.20.Cf;34.20.Gj;34.50.Lf
Keyword:transition state;quantum mechanics and molecular mechanics;solvent effect;potential of mean force
1. Introduction

Bimolecular nucleophilic substitution (SN2) reaction is a fundamental reaction in organic chemistry. In gas phase, the SN2 reactions have been extensively studied,[1–6] especially for the (X and Y are F, Cl, Br, and I) reactions.[7–12] In textbook, the standard SN2 reaction mechanism occurs with the attack of the Y− at the center carbon atom with a Walden inversion process. There is also a front-side attack reaction mechanism existing in these SN2 reactions.[13] Besides the Walden-inversion reaction mechanism, new reaction mechanisms have also been found, such as the rebound mechanism,[14] the roundabout mechanism,[15] the striping mechanism,[14] and the double-inversion mechanism.[16] Although the bimolecular nucleophilic SN2 reactions have been well studied in the gas phase, the studies of these reactions in solution are scarce.

In our previous studies of the CHmCl4 −m + OH−, CH3Br+OH−, CH3F+OH− in aqueous solution[17–21] with the multi-level quantum mechanics and molecular mechanics (ML-QM/MM) methods, the calculated reaction barrier heights agree well with the experimental results. The activation barriers increase significantly from gas phase to solution phase due to the presence of the aqueous solution, which indicates that the reactivity of the same reaction from gas-phase to solution is greatly reduced.

For the title reaction in gas phase, using the viral partitioning method, Bader et al. calculated an overall activation energy of 22.6 kcal/mol; the theoretical density functional theory (DFT) calculation with a 4-31G basis set gave the barrier height at 26.6 kcal/mol;[22] Ga’zquez[23] took into account of the difference of the hardness between the reactant and transition state and found the activation energy at 33.9 kcal/mol. In water, Marshall and Moelwyn–Hugh[24] measured the rate constants of the title reaction in aqueous solution and derived the activation energy of 23.0 ± 1 kcal/mol. We can see from the above studies that the activation barrier varies not much from the gas phase (26.6 kcal/mol) to aqueous solution (23.0 kcal/mol); somehow, the barrier height in gas phase is higher than the one measured in solution phase. This indicates that the water solution does not affect the reaction barrier height very much, even lowers the barrier height from gas phase to solution. This is very different from our previous calculations for other reactions[17–21] where the water solution raises the barrier height significantly from gas phase to solution.

Till now, no theoretical studies of the title reaction CH3F+CN− in aqueous solution have been done, therefore, in this article, employing the combined QM/MM method,[17–19] we carry out the first theoretical calculation of this reaction in water. In this study, the reactant solute region is treated with DFT[25] and electrostatic potential (ESP)[26] theories, and the solvent is treated with molecular mechanics. We want to study solvent effects on the reaction pathway, quantitatively determine the aqueous solution contribution to the transition state barrier, explore the detailed, atomic-level evolution of the reaction mechanism, obtain the potential of the mean force (PMF), and compare the activation barrier with the experimental result.[24]

2. Methods and computational procedures

In this article, the combined QM/MM method was used to study the reaction in aqueous solution. The reactant solute CH3F+CN− was described using the QM method, while the water environment was described using the MM method. The potential energy of the reaction system can be calculated as

Here and refer to the coordinates of the QM and MM parts, and ψ is the electron-wave function of the ground state of the solute. has the gas-phase expression and is the QM-region energy. The bond interaction between the nuclear solute and solvent, the van der Waals interactions, and the electrostatic interactions are included in the second term . The third term is the solvent energy contribution calculated using the molecular mechanics method.

Hence we can calculate the potential mean force, , between two neighboring points using the DFT/MM theory as

The last term is the solvent free energy difference calculated through explicit statistical sampling.

We first solvated the reactants CH3F and CN− into a 35.8 Å cubic box with 1532 SPC/E[27] solvent waters, treated as the MM region. The DFT/MPWB1K[25,28] theory with the aug-cc-pvDZ basis was used to describe the reactant QM part. The parameters for the van der Waals interactions were adopted from the standard Amber force.[29]

First, with a multi-region optimization procedure, the QM and MM regions were optimized. Hereafter, we fixed the QM part and implemented molecular dynamics simulation to equilibrate the solvent MM region for 40 picoseconds at 298 K. Second, we optimized the whole system again in water to obtain the initial reactant complex. Then, using the just-obtained reactant complex, we could seek the product complex according to the SN2 back-side attack mechanism. We also applied the multi-region optimization and equilibration procedures to obtain a converged initial product complex.

Next, the nudged elastic band (NEB)[30,31] method was employed to constructed the initial reaction path using the initial reactant and product complexes obtained above in water. We used the top stationary point from the initial NEB reaction path to search the real transition state in aqueous solution, following a numerical frequency calculation to confirm its saddle-point identity with one imaginary frequency.

We optimized the corresponding displacements of the imaginary frequency along the imaginary frequency mode. The NEB reaction pathway then was constructed anew with 10 points. For each of the NEB point, the water solution was equilibrated using molecular dynamics simulation for 40 picoseconds. Next, the NEB reaction path was optimized again. This equilibration-optimization procedure was repeated to reach a convergence. Finally, the PMF was calculated based on Eq. (2) on the converged NEB reaction pathway.

3. Results and discussion
3.1. Reactant complex

The comparison of the reactant complexes in gas phase[32] and aqueous solution is presented in Fig. 1. In gas phase, the C–C distance is 3.11 Å, while it is longer at 3.33 Å in aqueous solution. This is due to the shield effects resulting from the presence of the solvent between the nucleophile and the substrate, which reduces the interaction between the two. Another major difference is that the CCN angle is 180.0° in the gas phase, but non-linear in the aqueous solution. This feature is also due to the presence of the water molecules. This can be seen in that the CN− nucleophile accepts seven hydrogen bonds. Their average distance is 2.48 Å with the surrounding water molecules, and the uneven distribution of the water molecules around the nucleophile makes the CN− tilt up. This picture shows that the presence of the solvent changes the structure of the reactants in aqueous solution.

Fig. 1. (color online) Structures of the reactant complex for the reaction (a) in gas phase and (b) in aqueous solution. The indicated distances are in Å.
3.2. Transition state

The comparison of the transition states in gas phase[32] and in aqueous solution is shown in Fig. 2. The transition state [F…CH3…CN]− in water has one imaginary frequency at 425.7i cm−1. In gas phase, the C–C distance is 2.06 Å, while it is longer at 2.16 Å in aqueous solution; the CCN is 180.0° while it is 153.8° in aqueous solution. These two main differences are still due to the interactions between the solute and the surrounding water molecules as we discussed in the comparison of the reactant complexes above. The transition state has an almost planar CH3. However, the CN− gets closer to attack the C center of the CH3 group with the C–C distance at 2.16 Å, and we can see that the C–F distance is increased from 1.40 Å in reactant to 1.65 Å. In the meantime, seven water molecules form the solvation shell of the CN− nucleophile with an average distance at 2.45 Å.

Fig. 2. (color online) Structures of the transition for the reaction (a) in gas phase and (b) in aqueous solution. The indicated distances are in Å.
3.3. Product complex

The product complex is compared with the one in gas phase[32] in Fig. 3. In gas phase, the leaving group is totally detached from the C atom and forms a hydrogen bond with one of the H atoms in CH3 at 1.79 Å with the FHC angle of 172.2°. Meanwhile in aqueous solution, the leaving group behaves similarly as in gas phase and also forms a hydrogen bond with one of the H atoms in CH3, but the angle FHC is more bent at 148.5°. At this stage, the CN− nucleophile forms a new bond with the substrate with a C–C distance of 1.43 Å. There are still seven water molecules forming hydrogen bonds with CN at an average distance of 2.83 Å.

Fig. 3. (color online) Structures of the product complex for the reaction (a) in gas phase and (b) in aqueous solution. The indicated distances are in Å.

Note that the number of hydrogen bonds formed around CN/CN− is consistent with the results in previous studies.[33,34] Infrared and NMR study[35] indicates that the average hydrogen-bonded neighbors around the CN− include six to seven molecules, and molecular-dynamics calculations[34] show that the average number of hydrogen bonds formed with CN− in water is 6.3. Nonetheless, not all the hydrogen bonds formed with CN/CN− should be in its first solvation shell. According to the definition of weak hydrogen bond in the range of 2.0–3.0 Å,[35] the hydrogen bonds longer than 3.0 Å in the stationary points should be in the outer region of the CN/CN− solvation shell. Therefore, the CN/CN− solvation numbers of water molecules in the reactant, transition state, and product complex are six, five, and three, respectively. This is understandable as the CN− transfers its negative charge to the substrate, its electronegativity is decreased, thus, it loses water molecules in its first solvation shell.

3.4. Atomic-level reaction mechanism

In Fig. 4, we provide the atomic-level evolution of the geometries for the title reaction in water along the NEB reaction pathway. Ten snapshots are displayed along the reaction path in water. The C–C length decreases from 3.33 Å at reactant (Fig. 4(a)) to 2.16 Å at transition state (Fig. 4(e)); while in the meantime, the C–F distance increases from 1.40 Å to 1.65 Å. Moreover, the C–C distance further decreases from 2.16 Å at transition state to 1.43 Å at product state, and the C–F distance increases from 1.65 Å to 2.76 Å. On the other hand, at the beginning of the reaction process, the 3H atoms in CH3F face the CN− with the three HCC angles at 70.9°, 67.8°, 71.4° (Fig. 4(a)), respectively; then at the transition state (Fig. 4(e)), the three HCF angles become 92.3°, 104.6°, 79.4°, which indicates that the 3H atoms are almost perpendicular to the F–C bond; finally, at the product complex, they are inverted to face the leaving group F−. These above two aspects show a synchronized C–F bond-breaking and C–CN bond-forming mechanism (C–CN) accompanied with a Walden-inversion of the CH3 group; therefore, this reaction mechanism shows the back-side attach SN2 reaction mechanism for the CH3F+CN− reaction in aqueous solution.

Fig. 4. (color online) Structures of (a)–(j) the ten snapshots along the NEB reaction pathway for the reaction in aqueous solution. Panel (a) is the structure of the reactant state, panel (e) is the transition state, and panel (j) is the product state for the reaction. The indicated distances are in Å.
3.5. Charge distribution evolution

The electronic density of the solute was represented by the electrostatic potential (ESP) charges and the charge for each atom was fitted utilizing the Dunlap charge fitting method.[36] Figure 5 shows the evolution of the charge distribution along the NEB path. At point 1 (reactant complex), the CH3F is neutral with F atom having the negative charge of −0.36 and CH3 having the positive charge of +0.33, and the CN− has the whole negative charge of −1.0. At point 5 (transition state), the negative charge has been partially transferred to the F atom at −0.52 and the CH3 has position charge around +0.34. At point 6, the leaving group F has more negative charge −0.95 with the attacking group CNʼs charge reduced to −0.41. At point 10, the negative charge has been totally transferred to the leaving group F−, whereas the CN has negative charge −0.20, which keeps the product CH3CN neutral with charge +0.20. This process also shows a concerted charge transfer process from the attacking group to the leaving group.

Fig. 5. (color online) Charge distributions of the CH3, CN, and F along the NEB reaction pathway for the reaction .
3.6. Potential of mean force

The PMF obtained using the DFT/MM method and the solvent energy contribution (the last term in Eq. (2)) are shown in Fig. 6. We use the reactant complex (point 1) as the reference energy point. The transition state barrier height is 19.7 kcal/mol, agreeing with the experimental[24] one at 23.0 kcal/mol. The solvent energy contribution is −1.8 kcal/mol to the transition state, −30.6 kcal/mol to the product state. It indicates that the solvent energy lowers the barrier height. This situation is in agreement with previous studies that, in gas phase, the transition state height is found at 26.6 kcal/mol[22] and 33.9 kcal/mol,[23] while it is lower at 23 ± 1 kcal/mol[24] in aqueous solution.

Fig. 6. (color online) Comparison of the potential of mean force calculated at DFT/MM levels of theory and the solvation contribution using the reactant state (point 1) as a reference.

The solvent effects have two sources: the solvation energy and the polarization effect. In Fig. 7, comparing the gas-phase energies and the solute-internal energies, we can obtain the polarization effect to the PMF. The polarization effect results from the solute structure perturbed by the presence of the solvent. Using the ten geometries excluding the solute and solvent interactions, we can obtain the gas-phase potential energy; also by excluding the solution contribution, we can obtain the internal energy of the reaction system. By comparing these two lines, it is found that the polarization effect contributes 2.6 kcal/mol, 6.1 kcal/mol, and 7.2 kcal/mol to the reactant, transition state, and product state, respectively. The net polarization effect contributes 3.5 kcal/mol and 4.6 kcal/mol to the reaction free energy. As a sum combining these two effects, the solvent effects contribute 1.7 kcal/mol to the activation barrier and −26.0 kcal/mol to the product state.

Fig. 7. (color online) Comparison between gas-phase and solute internal energies along the NEB reaction pathway under the CCSD(T)/MM representation using the gas-phase energy of the reactant state as a reference.

We notice that, for the title reaction, the reaction barrier heights in the gas phase and in the aqueous solution are close; however, for other SN2 reactions we studied before, the reaction barriers in aqueous solution are much higher than those in gas phase. For example, for the CH3Cl+CN− reaction, the reaction barrier in gas phase is 11.5 kcal/mol,[37] while it is 20.3 kcal/mol in aqueous solution.[20] For the title reaction, the current study shows that the total solvent effects contribute only 1.7 kcal/mol to the barrier height; however, for the CH3Cl+CN−, the total solvent effects contribute a much bigger value 11.4 kcal/mol to the reaction barrier. Second, since F has a larger electronegativity than Cl does, the transition state conformation of [F…CH3…CN]− is much more compact than that of [Cl…CH3…CN]−:[20] the C–C and C–F distances are 2.16 Å and 1.65 Å for the former, while they are 2.44 Å and 2.24 Å for the latter. Thus, the solvent destabilizes the transition state of the former much less than the latter. These lead to the solvent effect having a smaller contribution to the title reaction than that to the CH3Cl + CN− reaction.

3.7. From gas phase to solution phase

Based on the early literature results of the free energies of solvation of CH3F (−0.22 kcal/mol),[38] CN− (−72.0 kcal/mol),[39] CH3CN and F− (−103.0 kcal/mol),[39] and the reaction energy in gas phase of −5.0 kcal/mol,[22] We can deduce the reaction free energy in aqueous solution at −39.7 kcal/mol (see Fig. 8), which agrees with our result of −43.5 kcal/mol. In Fig. 8, we draw the schematic plot of the stationary points (transition state at 26.6 kcal/mol[22]) of the reaction in gas phase and the calculated reaction free energy obtained based on the gas-phase data against our calculated PMF with the DFT/MM theory.

Fig. 8. (color online) Comparison between the potential energy profiles calculated along the NEB reaction pathway in gas phase (red) and in aqueous solution (blue) for the title reaction . The black numbers denote the experimental values and the green numbers denote the calculated values using gas-phase data.
4. Conclusion

The solvent effects and potential of mean force of the CH3F+CN− reaction in water are studied using the combined QM/MM approach. This is the first theoretical study of this reaction in aqueous solution. The interactions between the solute and solvent substantially affect the geometries of the stationary points. The PMF calculated with the DFT/MM theory has an activation barrier height at 19.7 kcal/mol, agreeing with the experimental one at 23.0 kcal/mol. This calculation shows that the solvent energy and the polarization effect combined contribute 1.7 kcal/mol to the activation barrier and −26.0 kcal/mol to the product state. The calculated reaction free energy −43.5 kcal/mol also agrees with the one estimated using gas-phase data.

Reference
[1] Wilbur J L Brauman J I 1991 J. Am. Chem. Soc. 113 9699
[2] Wolfe S Mitchell D J 1981 J. Am. Chem. Soc. 103 7694
[3] Hu W P Truhlar D G 1994 J. Am. Chem. Soc. 116 7797
[4] Tanaka K Mackay G I Payzant J D Bohme D K 1976 Can. J. Chem. 54 1643
[5] Graul S T Bowers M T 1991 J. Am. Chem. Soc. 113 9696
[6] Hase W L 1994 Science 266 998
[7] Viggiano A Morris R A Paschkewitz J S Paulson J F 1992 J. Am. Chem. Soc. 114 10477
[8] Parthiban S Oliveira G Martin J M 2001 J. Phys. Chem. 105 895
[9] Viggiano A Paschkewitz J S Morris R A Paulson J F 1991 J. Am. Chem. Soc. 113 9404
[10] Cyr D M Posey L A Bishea G A Han C Johnson M A 1991 J. Am. Chem. Soc. 113 9697
[11] Manikandan P Zhang J Hase W L 2012 J. Phys. Chem. 116 3061
[12] Sun R Davda C J Zhang J X Hase W L 2015 Phys. Chem. Chem. Phys. 17 2589
[13] Xie J Hase W L 2016 Science 352 32
[14] Mann D J Hase W L 1998 J. Phys. Chem. A 102 6208
[15] Mikosch J Trippel S Eichhorn C Otto R Lourderaj U Zhang J X Hase W L Weidemüller M Wester R 2008 Science 319 183
[16] Szabó I Czakó G 2015 Nat. Commun. 6 5972
[17] Yin H Wang D Y Valiev M 2011 J. Phys. Chem. 115 12047
[18] Wang T Yin H Wang D Y Valiev M 2012 J. Phys. Chem. 116 2371
[19] Xu Y Wang T Wang D Y 2012 J. Chem. Phys. 137 184501
[20] Xu Y Zhang J Wang D Y 2015 J. Chem. Phys. 142 244505
[21] Li C Niu M Liu P Li Y Wang D Y 2017 Chin. Phys. B 26 103401
[22] Shaik S Schlegel H B Wolfe S 1992 Theoretical Aspects of Physical Organic Chemistry New York Wiley p. 239
[23] Gázquez J L 1997 J. Phys. Chem. A 101 8967
[24] Marshall B W Moelwyn-Hughes E A 1965 J. Chem. Soc. 1312 7119
[25] Hohenberg P Kohn W 1964 Phys. Rev. 136 B864
[26] Valiev M Garrett B C Tsai M K Kowalski K Kathmann S M Schenter G K Dupuis M 2007 J. Chem. Phys. 127 051102
[27] Berendsen H J C Grigera J R Straatsma T P 1987 J. Phys. Chem. 91 6269
[28] Zhao Y Truhlar D G 2004 J. Phys. Chem. 108 6908
[29] Fox T Kollman P A 1998 J. Phys. Chem. 102 8070
[30] Henkelman G Uberuaga B P Jónsson H 2000 J. Chem. Phys. 113 9901
[31] Sheppard D Terrell R Henkelman G 2008 J. Chem. Phys. 128 134106
[32] Gonzales J M Sidney Cox R III Brown S T Allen W D Schaefer H F III 2001 J. Phys. Chem. 105 11327
[33] Eaton G Pena-Nuñez A S Symons M C R 1988 Faraday Discuss Chem. Soc. 85 237
[34] Ferrario M McDonald I R Symons M C R 1992 Mol. Phys. 77 617
[35] Desiraju G R Steiner T 2001 The Weak Hydrogen Bond: In Structural Chemistry and Biology Oxford Oxford University Press p. 10
[36] Dunlap B I Connolly J W D Sabin J R 1979 J. Chem. Phys. 71 3396
[37] Safi B Choho K Geerlings P 2001 J. Phys. Chem. 105 591
[38] Hine J Mookerjee P K 1975 J. Org. Chem. 40 292
[39] Pliego J R Riveros J M 2000 Chem. Phys. Lett. 332 597